Longitudinal Viscous Hydrodynamic Evolution for the Shattered Colour Glass 

Condensate 



Akihiko Monnai 1 '^ and Tetsufumi Hirano^Q 

1 Department of Physics, The University of Tokyo, Tokyo 113-0033, Japan 
(Dated: January 13, 2013) 

We investigate hydrodynamic evolution of the quark gluon plasma for the colour glass condensate 
type initial conditions. We solve full second-order viscous hydrodynamic equations in the longitu- 
dinal direction to find that non-boost invariant expansion leads to visible deformation on the initial 
rapidity distribution. The results indicate that hydrodynamic evolution with viscosity plays an 
important role in determining parameters for the initial distributions. 
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The heavy ion program at Large Hadron Collider 
(LHC) in European Organization for Nuclear Research 
(CERN) opens up new opportunities to explore the de- 
confined matter, the quark gluon plasma (QGP) [l|, in 
a wider temperature region. There would be also a 
good opportunity to investigate the colour glass con- 
densate (CGC) [3, i.e., a universal form of colliding 
hadrons/nuclei at very high energies, of which we had 
a glimpse at Relativistic Heavy Ion Collider (RHIC) in 
Brookhaven National Laboratory (BNL) 0]. 

Heavy ion reactions at high energies undergo several 
stages such as collisions of two nuclei, large entropy pro- 
duction just after the first contact, local thermalisation, 
hydrodynamic evolution and chemical/thermal freeze- 
outs of the system. Since a framework to describe the 
CGC is regarded as an effective theory of high energy 
hadrons/nuclei, it is often employed to describe the very 
first stage in high energy heavy ion collisions and to calcu- 
late multiplicity and/or rapidity distributions without as- 
suming secondary interactions [4-llQl]. Nevertheless, the 
data are remarkably described in this approach. On the 
other hand, relativistic hydrodynamic models are quite 
successful to describe space-time evolution of the QGP 
created in high energy heavy ion collision s, in pa rticular, 
anisotropy of transverse collective flow |ll| - [l4j . Initial 
states of the hydrodynamic evolution, however, are still 
uncertain so that quantitative conclusion about trans- 
port coefficients depends on initial modeling 15 1. This 
has been a longstanding issue which should be by all 
means resolved towards full understanding of the QGP. 
Both the initial gluon production from the CGC and the 
hydrodynamic evolution of the QGP are two distinct fea- 
tures of the whole reaction so that it is indispensable to 
unify these features [l|| and to dynamically model the 
reactions as a whole from the colliding two nuclei to final 
observablcs in a comprehensive fashion. 

A vast body of relativistic ideal and viscous hydrody- 
namic simulations have been performed so far to explore 
the bulk and transport properties of the QGP Boost 



invariant expansion in the longitudinal direction [1S|] is, 
however, often assumed in most of these simulations to 
reduce the numerical efforts even though boost invariant 
rapidity distributions have never been observed. Since 
particle production in low transverse momentum region 
is dominated by the small x modes in the nuclear wave 
function, where x is a momentum fraction of incident 
partons, the QGP production could be traced back to 
the initial partem density at small x inside the colliding 
nuclei before collisions. The gluon density at small x 
rapidly increases with decreasing x and is eventually sat- 
urated due to non-linear interactions among gluons. The 
non-boost invariant gluon production is predicted within 
the fc^-factorisation formula as a consequence of x depen- 
dence of saturation scale, Q 2 s {x) cx x~ x , where A ~ 0.3 
@, El controls rapidity and collision energy dependences 
of the saturation scale. 

In this Letter we describe hydrodynamic evolution of 
the hot matter in the beam direction with initial condi- 
tions from the shattered CGC. We focus especially on 
how hydrodynamic expansion with or without viscos- 
ity changes the flow-rapidity dependence of the entropy 
which could be identified with final rapidity distribution 
of hadrons. We neglect transverse hydrodynamic expan- 
sion in this study. Instead, we do not impose boost- 
invariance in the longitudinal expansion |20i422j . 

We solve the full second order constitutive equations 
23 1 generalised from the ones in the Israel- Stewart theory 
24j in the (l+l)-dimensional t-t} s coordinates. Here r 
and rj s are the proper time and the space-time rapidity, 
respectively, defined as t — r cosh rj s and z — t sinh r\ s . 
We neglect the baryon number current since we consider 
only gluon production in the /c^-factorisation formula as 
a initial condition. We choose the Landau frame where 
the energy dissipative current vanishes — 0. Then 
the constitutive equations are 
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Duhcm relation, the two linear terms are merged into one 
term at the first order as 



(2) 



where II is the bulk pressure and n the shear pressure de- 
fined with the shear stress tensor 7r Ml/ as tt 
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Note that, in the (l+l)-dimensional geometry, it is suf- 
ficient to treat 7r as the only independent component 
from the orthogonality and the traceless conditions. T 
is the temperature and Yf is the flow rapidity defined 
by the four-fluid velocity u M = (coshYf, 0, 0, sinhY/). r\ 
is the shear viscosity, £nn and Cn5e the bulk viscosities, 
Tn and T n the relaxation times and Xnn: Xnn> Xn>, X«i 
X%tt, xt-K and XnTi the second order transport coefficients. 
Here we have two "bulk viscosities" , £rm and £n<Se , as re- 
quired in the linear response theory. See also Ref. [23l j. 
The time- and the space-like derivatives in this geom- 
etry arc D = cosh(Y/ — f] s )d T + -sinh(Yf — rj s )d Vs and 
V = sinh(Y) — ij s )d T + -cosh(Y) — rj s )d ns . Hydrody- 
namic equations become complicated in non-boost invari- 
ant case Yf ^ rj s because derivatives with respect to the 
proper time and the space-time rapidity are mixed, which 
significantly increases the numerical difficulty compared 
with the constitutive equations in the transverse plane 
assuming the boost invariant flow. 

We solve the constitutive equations (p} and ((2]) to- 
gether with the energy-momentum conservation equa- 
tions in the piecewise parabolic method [25j which was 
employed in one of the most successful (3+1 )-dimensional 
ideal hydrodynamic calculations [l3[ . The numerical dif- 
ficulty raised by the mixing of the derivatives with respect 
to r and rj s is dealt with by taking iteration on the ex- 
pansion scalar 9 = VI/. We have checked the solutions 
converge typically in several steps. Numerical details will 
be discussed elsewhere j2r| . 

One needs the equation of state Pq = -Po(eo) and the 
transport coefficients, which contain microscopic infor- 
mation of the systems, to perform hydrodynamic simu- 
lations. We employ the latest (2+l)-flavor lattice QCD 
result [13] for the equation of state. On the other hand, 
we introduce some models for the transport coefficients 
since, to our knowledge, there is no single framework 
which gives all the transport coefficients appeared in this 
study. Here we use the conjectured minimum bound 
for the ratio of shear viscosity to the entropy density 
r]/s = 1/4-7T from Anti-de Sitter/conformal field theory 
(AdS/CFT) correspondence [28| just for the purpose of 
demonstration to see how visible the entropy distribu- 
tion changes during the longitudinal evolution. On the 
other hand, there are two bulk viscous coefficients £nn 
and CuSe when the first-order cross terms are properly 
kept. So far there have been few calculations on these 
coefficients. Hence we try to get an insight for these 
two coefficients from the </> 4 -theory in the non-equilibrium 
statistical operator method [29| and calculate the ratios 
Cnn/i] and QaSe/f] as 5T/6 and — 5T/2, respectively. Us- 
ing the energy-momentum conservation and the Gibbs- 
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where £ = (Cnn + c 2 s C,u5e)/T = |(| — c 2 )r] corresponds to 
the conventional bulk viscous coefficient. The relaxation 
times rn and tv and the other second order transport co- 
efficients Xnn> Xnn> Xnir, X™, X™,^lr and x^n are es- 
timated in kinetic theory as in Ref. [23| . The coefficients 
estimated with hadronic components up to A(1232) are 
connected to those with it, <i, s quarks and gluons as its 
components by hyperbolic factors around the (pseudo- 
critical temperature as x = ^(1— tanh )Xhad + ^(l+ 
tanh^o) Xudsg where T = 0.17 GcV and AT = 0.02 
GcV. The relaxation times are calculated likewise. We 
emphasize here that they arc trial parameters to inves- 
tigate qualitative response of the hot matter and that 
obtaining realistic transport coefficients is not the aim of 
this Letter. 

The initial conditions for hydrodynamic simulations 
are obtained from Monte Carlo version [30[ of kr- 
factorisation formula with unintcgratcd gluon distribu- 
tions parametrised by Kharzeev, Levin and Nardi (MC- 
KLN) [j| [3]} . In this model the saturation scale Q s for a 
nucleus A at a transverse coordinate x± is given as 
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where Ta(x±) is a thickness function which is obtained 
from randomly distributed hard source (nuclcons) ac- 
cording to the Woods-Saxon nuclear profile. We use 
the same parameter set, Q 2 S = 2 GeV 2 , Ta o = 1-53 
fm- 2 , A = 0.28 and x = 0.01 as in Ref. 0. The 
5% most central events in the Monte-Carlo calculations 
are employed for constructing the smooth initial condi- 
tions. The initial energy density distribution as a func- 
tion of space-time rapidity, eo{To,r] s ) is obtained from 
the transverse energy distribution dEx/dy over the over- 
lapping area of the nuclei >S ar ea by identifying momen- 
tum rapidity y with space-time rapidity rj s . Although 
the pre-thermalisation stage would be very complicated, 
we simply assume, during a very short period until the 
thermalisation time to = 1 fm/c, the energy density in 
a fluid element at tq is the same as locally deposited en- 
ergy dE /d 2 x±Todr/ s calculated in the CGC. Maximum 
(minimum) px cut is set to 3 (0.1) GeV/c to include 
the contribution from low p^ partons which is assumed 
to form thermally-equilibrated media. When the satura- 
tion scale is smaller than Aqcd, we simply assume gluons 
are not produced in fc^-factorisation formula. This pre- 
scription leads to reduction of the gluon multiplicity at 
mid-rapidity by 21.5% at RHIC and 39.9% at LHC. The 
temperatures at mid-rapidity at the initial time are 419 
MeV for RHIC and 490 MeV for LHC in these settings. 
We assume the boost-invariance only for the initial flow, 
i.e., Y/(to) = rj s , not for any other hydrodynamic initial 
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variables. It is more difficult to determine the initial con- 
ditions for the shear pressure 7r and the bulk pressure II 
because the precise dynamics before the local thcrmali- 
sation is not fully known. Here we set them to be vanish- 
ing since assuming finite dissipative currents at To corre- 
sponds to employing different initial energy-momentum 
tensors between ideal and viscous cases, which would 
make a comparison of the results and estimation of the 
viscous effects difficult. 

We investigate the entropy distribution as a function 
of the flow rapidity dS/dYf which roughly corresponds 
to the rapidity distribution of hadrons [33| . This can be 
interpreted as follows. The entropy density s carries the 
information of the number density in relativistic massless 
ideal gas limit since the ratio between the two quantities 
should be temperature independent as indicated by a di- 
mensional analysis. Momentum rapidity y, on the other 
hand, can be identified with the flow rapidity Yf in a fluid 
element on average. Thus, this is the closest quantity to 
the rapidity distribution we have from a pure hydrody- 
namic point of view, which do not include complicated 
freezeout processes nor any other model assumptions. 
Since the second order viscous correction to s M should not 
be large, we discuss here the lowest order modification. 
One can make a rough estimation of charged particle ra- 
pidity density by dN^ dm /dy w (2/3) x (1/3.6) x dS/dY f . 

The flow rapidity dependences of entropy are shown in 
Fig. [Q for Au+Au and Pb+Pb collisions at ^/snn = 200 
GeV and 2.76 TeV, respectively. The final times are r = 
30 fm/c and r = 50 fm/c for the Au+Au and the Pb+Pb 
cases, respectively. These are the typical times at which 
the temperature at mid-rapidity is sufficiently close to the 
(pseudo-)critical temperature T c ~ 0.17 GeV. Note that 
they are much longer than a conventional lifetime of the 
QGP in the heavy ion collisions because the transverse 
expansion neglected in the present study would accelerate 
the cooling process. 

One can see that the hydrodynamic evolution visibly 
modifies the initial distributions from the CGC model 
to more flattened ones at both RHIC and LHC energies. 
This behaviour can be seen only if boost- invariance is not 
assumed, i.e., a pressure gradient with respect to ij s ex- 
ists. Since entropy is produced in non-equilibrium hydro- 
dynamic evolution, the entropy distributions are larger 
than those of the ideal hydrodynamic systems in most 
of the flow rapidity region. Final entropy at mid-flow 
rapidity results from the interplay between the entropy 
production and outward flow which carries entropy to the 
forward rapidity region. The ideal hydrodynamic process 
always lowers the yield of the initial distribution at mid- 
flow rapidity due to the outward flow. The viscous hy- 
drodynamic one gives rather non-trivial results since the 
distribution is still lowered for the 200 GeV Au+Au colli- 
sions, but is enhanced for the 2.76 TeV Pb+Pb collisions. 
The difference is not due to the different choices of the 
final times, because, as we will see later, the shapes of the 
distributions do not change so much after t ~ 20 fm/c for 
both cases. The modification of the initial distribution 



is rather sensitive to the shape of the initial distribution. 
It should be noted that the corrections in the larger Yf 
region (Y f > 3 at RHIC and Y f > 6 at LHC) might be 
slightly overestimated because even though the distribu- 
tion shapes steady at a relatively early stage, the tem- 
perature in the forward region still hits T c earlier than 
the mid-rapidity regions do and there would be addi- 
tional hydrodynamic correction at the final times. This, 
of course, does not change the conclusion that the yield 
at the smaller rapidity regions are affected by hydrody- 
namic evolution. Considering that we use small shear 
and bulk viscous coefficients close to the conjectured min- 
imum boundaries [28|, [H| in the calculation, the results 
also suggest that effects of viscosity would be important. 
It should be noted that the numerical results here include 
only the longitudinal expansion and that the transverse 
expansion, which is missing in the present study, would 
be important in more quantitative discussion. 

The fact that the initial entropy distribution dS/dYf « 
dS/dy cx dN / dy could be deformed during the hydrody- 
namic stage is of particular importance since initial gluon 
rapidity distribution is often directly compared with 
the observed ch arged hadron data assuming the parton- 
hadron duality [4l-ll0l|. As we saw above, the hydrody- 
namic expansion beyond the boost invariant flow could 
make a change of initial rapidity distribution. Thus, the 
A parameter in some CGC-inspired models is subject to 
correction due to hydrodynamic effects if one wants to 
constrain it through the rapidity distribution. Even if 
one satisfactorily reproduced the rapidity distribution at 
RHIC with a set of CGC parameters assuming parton- 
hadron duality, one could fail to reproduce LHC data 
with the same parameter set due to lack of hydrodynamic 
corrections. Suppose the result from the viscous hydro- 
dynamic model in Fig. [T] (a) would reproduce hadron 
rapidity distribution at RHIC, then the actual A which 
includes the hydrodynamic expanding effects would be 
larger than an apparent A which can be obtained by fit- 
ting this rapidity distribution without considering any 
secondary scatterings, as the distribution tends to be- 
come steeper as increasing A. If one assumes that the 
hydrodynamic effect is smaller at LHC as indicated in 
Fig. Q] (b), this could be one of the possible interpreta- 
tions for the fact that most of CGC models turned out to 
underpredict the multiplicity in the Pb+Pb collisions at 
LHC [35[ . Since the entropy production due to viscosity 
is non-trivial as increasing collision energy, the energy 
dependence of multiplicity even at mid-rapidity from the 
CGC is also subject to hydrodynamic correction. 

To further quantify the deviation from the boost invari- 
ance (Yf = r) s ), we also estimate the difference between 
the flow rapidity and the space-time rapidity, Yf — rj s . 
Results for Au+Au collisions at ^/snn = 200 GeV and 
Pb+Pb collisions at ^/snn = 2.76 TeV are shown up to 
the beam rapidities in Fig. [2j One sees in all cases that 
the deviations from the boost-invariant flow are positive 
and become large towards larger space-time rapidity due 
to the acceleration by the pressure gradient in the beam 
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FIG. 1. (Colour online) The initial dS/dYj distributions at t = 1 fm/c from the colour glass condensate (solid line) and 
the final distributions after the ideal hydrodynamic (dotted line) and the shear and bulk viscous hydrodynamic (dashed line) 
evolution for (a) Au+Au collisions at RHIC (r = 30 fm/c) and (b) Pb+Pb collisions at LHC (r = 50 fm/c). 
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(b) initial flow, Pb+Pb 2.76TeV 
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FIG. 2. (Colour online) The deviation of the flow rapidity from the space-time rapidity Yf — r/ s at the initial time r = 1 fm/c 
(solid line) and after the ideal hydrodynamic (dotted line) and the shear and bulk viscous hydrodynamic (dashed line) evolution 
for (a) Au+Au collisions at RHIC (r = 30 fm/c) and (b) Pb+Pb collisions at LHC (r = 50 fm/c). 



direction. Both Au+Au case at RHIC and Pb+Pb case 
at LHC exhibit the same trend, while the latter is slightly 
moderate near mid-rapidity possibly because initial en- 
ergy density profile is less steep in the case. The devia- 
tions for the viscous cases are slightly smaller than those 
for the ideal cases for the Yf regions in which the fluids 
are relatively hot, because the pressure gradients are ef- 
fectively reduced in the longitudinal direction by the bulk 
and the shear pressures at early times in the space-time 
evolution. The situation, however, is different for the 
later times because the entropy generated in the viscous 
systems enhances the pressure gradients, while the cor- 
rections from the shear and the bulk pressures themselves 
are already small. Due to the counter contributions at 



the late stage, the overall differences between the ideal 
and the viscous flow deviations are small at those proper 
times. 

Finally, we demonstrate the time evolution of dS/dYf 
and Yf — r] s in Pb+Pb collisions at y/s^N = 2.76 TeV. 
In Fig. [3] (a), the proper time dependence of the entropy 
distribution per flow rapidity is shown at the time r = 1 , 
5, 20 and 50 fm/c. As mentioned earlier, the entropy dis- 
tribution in the hydrodynamic evolution does not change 
much its shape after r ~ 20 fm/c. The yield around 
Yf = is almost constant throughout the time evolution. 
It is due to a rather accidental cancellation between the 
entropy production and the expansion by the outward 
entropy flux in the parameter settings. This means we 
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FIG. 3. (Colour online) (a) The deformation of the initial entropy distribution per flow rapidity and (b) the deviation of the 
flow rapidity from the space-time rapidity Yf — rj a at initial time r = 1 fm/c (solid line), r = 5 fm/c (dotted line), r = 20 fm/c 
(dash-dotted line) and r = 50 fm/c (dashed line) in viscous hydrodynamic evolution in Pb+Pb collisions at the LHC energy. 



have a monotonous decrease at mid-flow rapidity in the 
case of an ideal hydrodynamic calculation. In Fig. [3] (b) , 
one sees the dynamical evolution of the deviation from 
the boost-invariant flow. At r = 5 fm/c, there are an rise 
and a dip in the flow acceleration near Yf ~ 6-7 and 7-8, 
respectively, because the effective pressure Po+H — n can 
become very small in the large Yf region when the abso- 
lute values of II and tt are still large. The sudden decrease 
in pressure leads to an enhancement in its gradient fol- 
lowed by a suppression. Note that the effect quickly dis- 
appears as the dissipative currents rapidly approach van- 
ishing along with the time evolution. Eventually the flow 
rapidity distribution evolves into the one closer to the 
ideal hydrodynamic distribution we have seen in Fig. [2] 
It is worth- mentioning that, unlike dS/dYf, the flow ra- 
pidity profile changes in the time evolution after t = 20 
fm/c for the current parameter sets. 

To summarize, we developed a (l+l)-dimensional 
second-order viscous hydrodynamic model with both 
shear and bulk viscosity to see the QGP dynamics in 
the longitudinal direction. There is no boost-invariance 
at both RHIC and LHC energies in the CGC initial 
conditions, which causes the deformation of the entropy 
per flow rapidity. This indicates that the shapes of the 
(pseudo-) rapidity distributions of hadrons observed in ex- 
periments would reflect the initial gluon rapidity distri- 



butions only indirectly due to the hydrodynamic evolu- 
tion. This also motivates ones to correct the parameters, 
in particular, A which controls the rapidity dependence of 
entropy production, in the initial conditions. While the 
precise determination of the parameters should be left to 
a (3+l)-dimcnsional viscous hydrodynamic calculation, 
our current parameter settings suggest that a conven- 
tional A which is determined without the hydrodynamic 
effect would be smaller than the true A at RHIC but is 
not so much different at LHC around the mid-rapidity 
region. This could play an important role in explain- 
ing the gap between the CGC predictions of multiplicity 
based on RHIC data and the latest experimental data 
at LHC [35[. In future, more realistic wave functions 
under running-coupling quantum evolution will be em- 
ployed (36[. We will discuss parameter dependences and 
numerical aspects in detail in near future [26| . 
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